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We derive a closed equation ol motion for the current density of an inhomogeneous quantum many- 
body system under the assumption that the time-dependent wave function can be described as a 
geometric deformation of the ground-state wave function. By describing the many-body system in 
terms of a single collective field we provide an alternative to traditional approaches, which emphasize 
one-particle orbitals. We refer to our approach as continuum mechanics for quantum many-body 
systems. In the linear response regime, the equation of motion for the displacement field becomes 
a linear fourth-order integro-differential equation, whose only inputs are the one-particle density 
matrix and the pair correlation function of the ground-state. The complexity of this equation remains 
essentially unchanged as the number of particles increases. We show that our equation of motion is 
a hermitian eigenvalue problem, which admits a complete set of orthonormal eigenfunctions under a 
scalar product that involves the ground-state density. Further, we show that the excitation energies 
derived from this approach satisfy a sum rule which guarantees the exactness of the integrated 
spectral strength. Our formulation becomes exact for systems consisting of a single particle, and for 
any many-body system in the high-frequency limit. The theory is illustrated by explicit calculations 
for simple one- and two-particle systems. 



I. INTRODUCTION 

The dynamics of quantum many-particle systems 
poses a major challenge to computational physicists and 
chemists. In the study of ground-state properties one 
can rely on a variational principle, which enables a vari- 
ety of powerful statistical methods (in addition to exact 
diagonalization) such as the quantum variational Monte 
Carlo method and the diffusion Monte Carlo method. 1 
In time-dependent situations, the absence of a practi- 
cal variational principle has greatly hindered the devel- 
opment of equally powerful methods. Yet it is hard to 
overestimate the importance of developing effective tech- 
niques to tackle the quantum dynamical problem. Such 
a technique could allow, for example, to follow in real 
time the evolution of chemical reactions, ionization and 
collision processes. 

One of the most successful computational meth- 
ods developed to date is the time-dependent density 
functional theory (TDDFT), or its more recent ver- 
sion - time-dependent current density functional theory 
(TDCDFT). 2 In this approach, the interacting electronic 



system is treated as a noninteracting electronic system 
subjected to an effective scalar potential (a vector poten- 
tial in TDCDFT) which is self-consistently determined 
by the electronic density (or by the current density). 3 ^ 
Thus, one avoids the formidable problem of solving the 
time-dependent Schrodinger equation for the many-body 
wave function. Even this simplified problem, however, 
is quite complex, since it involves the determination of 
N time-dependent single particle orbitals - one for each 
particle. Furthermore, there are features such as multi- 
particle excitations 5 and dispersion forces 6 that are very 
difficult to treat within the conventional approximation 
schemes. 

An alternative approach, which actually dates back 
to the early days of the quantum theory, attempts 
to calculate the collective variables of interest, density 
and current, without appealing to the underlying wave 
function^— This approach we call "quantum continuum 
mechanics" (QCM), because in analogy with classical 
theories of continuous media (elasticity and hydrodynam- 
ics) , it attempts to describe the quantum many-body sys- 
tem without explicit reference to the individual particles 
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of which the system is constituted^ 

That such a description is possible is guaranteed by the 
very same theorems that lie at the foundation of TDDFT 
and TDCDFT. 11 ' 12 Indeed, consider a system of particles 
of mass m described by the time-dependent hamiltonian 

H(t) = H + J drn(r)^(r,i) (1) 

where 

H =f + W + V (2) 

is the sum of kinetic energy (T), interaction potential 
energy (W), and the energy associated with an external 
static potential Vo(r), 

% = J drVb(r)n(r), (3) 

where h(r) is the particle density operator. Vi(r,t) is an 
external time-dependent potential. 

The exact Heisenberg equations of motion for the den- 
sity and the current density operators, averaged over the 
quantum state, lead to equations of motion for the aver- 
age particle density n(r, t) and the average particle cur- 
rent density j(r, t): 

d t n(r,t) = -a M j M (r,i) (4) 

and 

md tJf t(r,t)= - n{r,t)d f ,[V Q (r) + V 1 (r,t)} 

- d v P„ u (r,t), (5) 

where d t denotes the partial derivative with respect to 
time and d u is a short-hand for the derivative with re- 
spect to the cartesian component v of the position vector 
r. Here and in the following we adopt the convention 
that repeated indices are summed over. These equations 
simply express the local conservation of particle number 
(Eq. flU) and momentum (Eq. (jSJ). The key quantity on 
the right hand side of Eq. ([5]) is the stress tensor P M i/(r, t) 
- a symmetric tensor which will be defined in the next 
section as the expectation value of a hermitian operator, 
and whose divergence with respect to one of the indices 
yields the force density arising from internal quantum- 
kinetic and interaction effects. 

Now the Runge- Gross theorem of TDDFT guarantees 
that the stress tensor, like every other observable of the 
system, is a functional of the current density and of the 
initial quantum state. Thus, Eq. ([3]) is in principle a 
closed equation of motion for j - the only missing piece 
of information being the explicit expression for P^ v in 
terms of the current density. 

In recent years much effort has been devoted to con- 
structing an approximate QCMi^r— and several applica- 
tions have appeared in the literature (see Ref. [22| for some 
representative examples). All approximation schemes so 
far have been based on the local density approximation 



and generalizations thereof. The objective of this pa- 
per is to derive and discuss an approximate expression 
for P M „(r, t), and, more importantly, for the associated 
force density —d v P^ u (r,t) : as functionals of the current 
density and the initial state. We will do this in the lin- 
ear response regime, i.e. for systems that start from the 
ground-state of the static hamiltonian Hq and perform 
small-amplitude oscillations about it. The external po- 
tential Vi(r, t) will be treated as a small perturbation. 
In this regime the equations of motion (Q} and (J5J) are 
conveniently expressed in terms of the displacement field 
u(r, t), defined by the relation 

j(r,i) = n (r)d t u(r,t) , (6) 

where no(r) is the ground-state density. It is also conve- 
nient to write the density and the stress tensor as sums 
of a large ground-state component and a small time- 
dependent part, in the following manner 

n(r,t) = n (r) + m(r,i) , 
P^(r,t) = P w o(r)+P^,i(r,t), (7) 

where the equilibrium components, marked by the sub- 
script 0, satisfy the equilibrium condition 

no(r)^Vb(r) + d u P^ fi (v) = . (8) 

Then the two equations (|4| and ([5]) take the form 

m(r, t) = -9,i[no(r)« M (r, t)} , (9) 

and 

mra (r)d t 2 u M (r,i) = -n (r)9 M Fi(r, t) 
-ni(r,t)a^o(r)-ft,P wX (r,t). (10) 

Our task is to find an expression for the force density 
5i/P /J1 /,i(r, t) as a linear functional of u(r, t). If this can 
be achieved, then the excitation energies of the system 
will be obtained from the frequencies of the time-periodic 
solutions of Eq. (ITU|) in the absence of external field (i.e., 
with Vi = 0). 

It is easy to see that the spatial dependence of these so- 
lutions will be proportional to the matrix element of the 
current density operator between the ground-state and 
the excited state in question. This is because, in a many- 
body system with stationary states \ipo}, fyi), ■■■■> \ipn)i ■■■ 
(\ipo) is the ground-state), and corresponding energies 
Eq, E\, E n ..., the n-th linear excitation is described 
by the time-dependent state 

\i) Q )e- iEat + e\i) n )e- iE "\ (11) 

where e is an arbitrarily small "mixing parameter" . The 
expectation value of the current density operator in this 
state is 

j(r,t) = e(V ; o|j(r)^n)e- <( ^- Eo)t + c.c. (12) 
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Thus, in principle, almost all the excitation energies 
(E n — E ) of the system can be obtained by Fourier- 
analyzing the displacement field - the only exception 
being those excitations that are not connected to the 
ground-state by a finite matrix element of the current- 
density operator. 

In this paper we will introduce an approximate expres- 
sion for the force density 

F M ,i(r,t) = -ni(r,t)S„%(r)-ft,P w i(r,t), (13) 

which appears on the right hand side of Eq. (ITU1) , as a lin- 
ear functional of u(r, t). The expression will be presented 
in terms of the functional 

B[u] = Wu]|ioWu]), (14) 

which is the energy of the distorted ground-state |-0o[ u ]) ; 
obtained from the undistorted ground-state |^o) by vir- 
tually displacing the volume element located at r to a 
new position r + u(r, t). More precisely, we will show 
that the equation of motion for u takes the form: 

mn (r)9 t 2 u(r,t) = -n (r) Wi(r, t) - ^| , 

ou(r, t) 

(15) 

where ^[u] is the second order term in the expansion 
of E[u] in powers of u. The functional ^[u] has an 
exact expression in terms of the one-particle density ma- 
trix and the pair correlation function of the ground-state, 
which is a major simplification, since ground-state prop- 
erties, unlike time-dependent properties, are accessible 
to computation by a variety of numerical and analytical 
methods. 

Furthermore, we will show that the kinetic part of the 
force density functional <5-E 2 [u]/<5u is local, in the sense 
that it depends only on a finite number of spatial deriva- 
tives (up to the fourth) of the displacement field at a 
given position. Thus, our equation of motion reduces to 
a fourth-order differential equation for u when interac- 
tion effects are neglected. The inclusion of interaction 
effects leads to the appearance of nonlocal contributions 
to the energy, and the equation of motion becomes a 
fourth-order integro- differential equation for the displace- 
ment field. However, the complexity of this equation re- 
mains essentially unchanged as the number of particles 
increases. 

Our equation of motion has two especially appealing 
features: (i) it is exact for one-electron systems at all fre- 
quencies and (ii) it can be physically justified for generic 
many-electron systems at high frequency or, more gener- 
ally, at all frequencies for which a collective description 
of the motion is plausible. Thus the range of frequencies 
for which our approximation makes sense is expected to 
be wider in strongly correlated systems than in weakly 
correlated ones. 

We discuss several qualitative features of our equation 
(uniform electron gas limit, harmonic potential theorem) 



and present its solution in simple one- and two-electron 
models, where the results can be checked against exact 
calculations. The results are encouraging. Although we 
are not able to resolve all the different excitation ener- 
gies of the models under study, we find that groups of ex- 
citation characterized by similar displacement fields are 
represented by a single mode of an average frequency, 
in such a way that the spectral strength of this mode 
equals the sum of the spectral strengths of all the excita- 
tions in the group. In this sense our approximation can 
be viewed as a (considerable) refinement and extension of 
the traditional single-mode approximation for the homo- 
geneous electron gas to strongly inhomogeneous quantum 
systems. In spite of the somewhat limited range of valid- 
ity of the present treatment (the linear response regime), 
we feel that this is an important first step in a direc- 
tion that might eventually lead to the construction of 
useful force density functionals for far-from-equilibrium 
processes. 



This paper is organized as follows. In Section II we 
present a complete derivation of the linearized equation of 
motion for the displacement field. We begin by deriving 
a formally exact expression for the force density (Section 

II A), on which we perform the "elastic approximation" 
(Section II B). The expression for the force density in the 
elastic approximation is worked out in sections II C (ki- 
netic part) and II D (potential part). A simplified form 
of the equation of motion, valid for one-dimensional sys- 
tems, is presented in Section II E. Appendixes A through 
C provide supporting material for this part. In Section 

III we discuss the relation between quantum continuum 
mechanics and time-dependent current density functional 
theory. In section IV we show how the linear equation 
of motion derived in Section II leads to an eigenvalue 
problem for the excitation energies. In Section IV A 
we demonstrate the hermiticity of this eigenvalue prob- 
lem and the positive-definiteness of the eigenvalues. In 
Section IV B we connect the eigenvalue problem to the 
high-frequency limit of the linear response theory. In 
Section IV C we prove that the first moment of the cur- 
rent excitation spectrum obtained from the solution of 
our eigenvalue problem is exact. Appendixes D and E 
contain supporting material for this part. In Section V 
we present a few simple applications of our theory for the 
excitations of (i) a homogeneous electron gas (Section V 
A) (ii) the linear harmonic oscillator and the hydrogen 
atom (Section V B), and (hi) a system of two-electrons 
in a one-dimensional parabolic potential interacting via 
a soft Coulomb potential. The analytic solution of the 
last model in the strong correlation regime is featured in 
Appendix F. Finally, Section VI contains our summary 
and a few speculations about future applications of the 
theory. 
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II. LINEARIZED EQUATION OF MOTION 
A. Derivation of the force density 

In this section we undertake the construction of an ap- 
proximate expression for the force density, Eq. (|13[) . as a 
linear functional of u. The stress tensor P^(r,t), whose 
divergence determines the force density, is defined as the 
expectation value of the stress tensor operator P^ u (r) in 
the evolving quantum state \ip{t)): 



p^t) = m)\h^)\m) 



(16) 



An exact and unambiguous expression for the operator 
Pfiv (r) hi an arbitrary system of coordinates is obtained 
by considering the universal many-body Hamiltonian 



H u =T + W 



(17) 



(external potential not included) in the presence of a 
"metric tensor field" g M „(r). As is well known^, the 
metric tensor g^ v (r) allows us to express the length ds of 
an infinitesimal displacement from r to r + dr in terms 
the corresponding increments of the coordinates dr^ : 



ds z 



g^vi^dr^dry . 



(18) 



In ordinary Euclidean space and in cartesian coordi- 
nates gfiv{?) = 5fj,vi independent of position. In gen- 
eral, however, a non-Euclidean space is characterized by a 
position-dependent, symmetric g^ v {r). A non-Euclidean 
metric can also be generated by a change of coordinates 
in an Euclidean space, as we will see shortly. As an im- 
portant technical point we also introduce the tensor g^ v 
as the inverse of and we define g as the determinant 
of g^v (so g 1 is the determinant of g^ u ). 

The hamiltonian H u undergoes the following changes 
in the presence of a non-trivial metrics. First, the lapla- 
cian operator d^d^ for the kinetic energy is replaced by 
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(19) 



Second, the Euclidean distance between two points 
(which controls the interaction energy) is replaced by the 
non-Euclidean length of the shortest path (geodesic) con- 
necting the two points. We denote by i? u [g] the Hamil- 
tonian in the presence of the metric field g^ v . Then 
the stress tensor operator is defined as the first varia- 
tion of H u [g] under an instantaneous virtual variation of 
the metric tensor (^(r), i.e., 



2 6H U [ S ] 



(20) 



The first-order change in the Hamiltonian due to a change 
Sg^n in the metric tensor is given by 



W, g //,-g] ! / dx^p-P^td&gTt?) . (211 



Notice that the stress tensor operator defined in this man- 
ner is itself a functional of the metrics. This definition 
is completely analogous to the standard definition of the 
current density operator as the derivative of the Hamil- 
tonian with respect to a vector potential. An explicit 
expression for P M „ in Euclidean metrics is reported for 
completeness in Appendix A (see also Refs [lH 
andl25l). We note that the definition of the quantum me- 
chanical stress tensor via the variational derivative with 
respect to the metric tensor has been also employed in 
Ref. H 

We will now focus on the calculation of P^i ~ the 
correction to P^ of first order in u. In order to ex- 
press Fylfot) and its divergence as functionals of the 
displacement field we resort to Tokatly's recent formu- 
lation of quantum dynamics in the co-moving reference 
frame.— ~— The co-moving frame is an accelerated refer- 
ence frame which, at each point and each time, moves 
with the velocity of the volume element of the fluid at 
that point and that time, so that the density is constant 
and equal to the ground-state density, while the current 
density is zero. The time-dependent transformation from 
the laboratory frame (coordinates r) to the co-moving 
frame (coordinates £) is defined by the solution of the 
equation 



d t r{t) =v(r,t), 



r(0)=C: 



(22) 
In the 



where v(r, t) = jjft^k is the velocity field, 
linear response regime the velocity is approximated as 
j(r, i)/no(r), where no(r) is the ground-state density. In 
this regime we can write 



r(t)=C + u(€,t) 



(23) 



where u(£,i) is the (small) displacement of a fluid ele- 
ment for its initial position Expressing ds 2 = dr ■ dr 
in terms of the new coordinates £ and making use of 
Eq. (fT8|) we see that the metric tensor in the co-moving 
frame is given by 



ft A - dra dra 



(24) 



From Eq. (|2"3")l . to first order in the displacement field, we 
immediately get 



and 



where 



or = <v - 2 Uflu 



(25) 
(26) 

(27) 



is the strain tensor. Also to first order in u the determi- 
nant of the metric tensor is easily seen to be 



g ~ 1 + 2V • u , 



(28) 
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so that, for example, g 1 I 2 ~ 1 — V • u. In view of these 
relations we will, in the rest of this paper, replace H u [g] 
by H u [u] , with the understanding that u completely de- 
termines the metrics. We also notice that, by virtue of 
Eq. (j2"6]l , we have 



6H u [g\ _ 1 6H u [u] 
StTir) 2*u pi/ (r) 



(29) 



The main reason for introducing the co-moving refer- 
ence frame is that in this frame we can make a simple 
approximation, which enormously simplifies the task of 
linearizing the stress tensor. This will be discussed in the 
next section. For the time being we proceed in a formally 
exact manner. To begin with, we observe that the gen- 
eral relation between the stress tensor in the lab frame 
and that in the co-moving frame is 



P MV (r,t) = (d^ a )(d^p)P a ^(r,t),t) , (30) 



where 



iV(r,i) = - 



-(m\P^\m), (si) 



where \ij)(t)) is the quantum state in the co-moving 
framed 

After expanding the stress tensor in the co-moving 
frame to first order in the displacement field, 

P^id-t) = +iW(€,t) . (32) 

it is easy to see that in the lab frame we have 

Ppv,l = P^lvA - U • VP^O - (d,j.Ua)Pav,0 ~ {d v U a ) P^afi ■ 

(33) 

From this we get 

9 v Pfj,v,l — dv[Pfj,v,l — (d^Ua)P av fl — (duU^P^afl] 

- a„(u-VP F , ), (34) 
and, after some algebra, 

nid^Vo + dvP^vA = n u ■ VO^Vq 

+ d v [Pnv,l + (V • u)P M „,o - (dfiV,a)Pau,0 ~ ^U av P^, a fi] 

(35) 

where we have made use of the equilibrium condition ([8jl 
and the definition (j2"?| of the strain tensor. 

It is convenient at this point to introduce the first-order 
stress force density 

= -dv[P^u,i+{^-^)PiJ,u,a-{d IJ ,u a )P aV fi-2u au P lla fi} 

(36) 

so that the equation of motion (flUt takes the form 



mdfu^v, t) + u ■ Vd,V = _ dVl ( T , t) . 

n (r) 



Finally, it is possible to prove (see Appendix B) that 
the first order force density is exactly given by the ex- 
pression 



J>,i(r,t) = -Wt)|^M|#)) 



(38) 



where 



8H U 



is its functional derivative calculated with 



respect to a virtual (time-independent) variation of the 

displacement field. The vertical bar mandates that 

l 

we keep only the first-order in u part of the bracketed 
expression. 

With the help of this identity we see that the equation 
of motion for the displacement field takes the form 



mdfu^T^) + u- Vdf.Vo 

1 



n (r) 



(39) 



The same result could have been derived almost imme- 
diately by using the more sophisticated machinery of the 
generally covariant Lagrangian formalism introduced in 
Ref. In fact, Eq. (|39[) is simply a linerized version of 
the equation of motion for an infinitesimal fluid element, 
Eq. (39) of Ref. US 

As a reality check, let us ask ourselves whether the sys- 
tem can support excitations in which the displacement 
field is uniform in space: u(r, t) = u(t). Clearly in this 
case the strain vanishes and there is no change in met- 
rics, so SH u [u]/Su is null. As a result, after setting the 
external field V\ = we get the equation 



(40) 



which has the solution u^{t) cx cos{ujt + <f>) if and only if 
the potential is of the harmonic form Vb(r) = ^muj 2 r 2 . 
This is just a statement of the harmonic potential theo- 
rem*^ according to which a many-body system in a har- 
monic potential performs a rigid simple harmonic motion 
with frequency uj imposed by the curvature of the har- 
monic potential. We have now shown that the harmonic 
potential is the only potential with this property. 



B. The elastic approximation 



(37) 



Equation (|39|) is formally exact, but it still contains 
the time-dependent state \ip(t)), which of course is not 
known. In spite of the simple behavior of the density 
(constant) and the current density (null), the evolution 
of the many-body wave function in the co-moving frame 
is far from trivial. Nevertheless, a simple and physically 
appealing approximation suggests itself. Namely, we as- 
sume that the wave function -0 is time-independent (just 
as the density) , and coincides with the ground-state wave 
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function of the laboratory frame (i/'o) evaluated at the co- 
ordinates £ of the co-moving frame: 



(41) 



The physical idea behind this approximation is that 
the time evolution of the wave function in the laboratory 
frame can be approximated as a continuously evolving 
elastic deformation of the ground-state wave function. 
Such a deformation affects all the particles simultane- 
ously and instantaneously. The burden of describing the 
time-evolution of the system is entirely placed on the 
time-dependent geometry (i.e. the time-dependent re- 
lation between £ and r), while the wave function itself 
remains independent of time. 

What is lost in this approximation is the fact that in 
the actual time evolution the system will undergo internal 
relaxation in order to optimize the correlations between 
the particles. In other words, the probability of finding 
the particles in a certain configuration ri, ...rjv at time t 
is not strictly determined by the probability that those 
particles were initially in the configuration £1, £jv from 
which r!,...rjv evolve according to Eqs. (f2"2"j) and ([23]). 
However, our approximation should always be valid at 
sufficiently high frequency, i.e., when the evolution of the 
geometry is very fast on the scale of the characteristic 
response times of the system. 

The equation of motion resulting from the elastic ap- 
proximation is also strictly valid (and therefore, not an 
approximation at all) for any one-particle system, be- 
cause in this case the wave function is completely deter- 
mined by the displacement field and there is no room 
for internal relaxation. Finally, our equation of motion 
is also strictly valid for non-interacting Bose systems in 
the ground-state (since these systems behave like a single 
particle), and for non-interacting Fermi systems consist- 
ing of at most two particles of opposite spins in the same 
orbital (since these behave like non-interacting Bosons). 
In all other cases - including the apparently simple case 
of a non-interacting many-fermion system - the appro- 
priateness of the elastic approximation must be assessed 
a posteriori and may depend on the objective of the cal- 
culation. In general, we can only say that the elastic 
approximation is expected to work better for collective 
(many-particle) excitations than for single particle ex- 
citations, and better for strongly correlated many-body 
systems (which exhibit bosonic behavior) than for weakly 
correlated systems. 

It is important to appreciate the profound difference 
that exists between the present approximation and an- 
other common approximation which also entails an in- 
stantaneous response to a time-dependent field: the adia- 
batic approximation. In the adiabatic approximation one 
assumes that the system remains in the instantaneous 
ground-state of the hamiltonian H(t) - an assumption 
that is justified only if the time evolution is slow on the 
scale of the characteristic response time of the system. 
This is exactly the opposite of the regime of validity of 
the present approximation. The geometrically distorted 



wave function ip is not at all close to the instantaneous 
ground-state of H (t) . Rather, it is the ground-state of the 
"deformed hamiltonian" i?o[u] which is obtained from 
the initial-time Hamiltonian Hq by a coordinate trans- 
formation - indeed an elastic deformation. 

As anticipated in the foregoing discussion the elastic 
approximation paves the way for a relatively simple cal- 
culation of the complicated expression that appears in 
the second line of Eq. (|3T)1) . Namely, thanks to the fact 
that ip = tpo is independent of the displacement field we 
can take the functional derivative of Eq. after taking 
the average and we arrive at 



•7>,i(r,t) 



SEJu] 



where 

E u [u] = (^ \H u [u]\^ ) . 
We further observe that 

1 SV [u] 



u-Vd^V (r) = 



n (r) (5M M (r) 



where 



V [u] = J drF (r + u(r))no(r). 



(42) 



(43) 



(44) 



(45) 



is the external potential energy of the distorted ground- 
state. 

Putting all together we arrive at the elegant result: 



1 SE 2 [u] 
n (r) <Su M (r) 



-3 M Vi(r,t). 



(46) 



where E2 [u] is the second-order term in the expansion of 
the total energy 



E[u] = E u [u] + V [u} 



(47) 



of the distorted ground state. Equivalently, E[u] can be 
obtained as the expectation value of the original hamil- 
tonian H in the distorted ground state 



N 



where & = r^ — u(r-j) and the last factor on the right hand 
side is for normalization. This proves that E[u] — E[0] is 
a positive definite quantity since E[0] is the ground-state 
energy of Hq while E[u] is the expectation value of Hq 
in a state that is not the ground-state. 
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C. Calculation of 8E'2[u]/Su — kinetic part only a finite number of derivatives (up to the fourth) of 

the displacement field. For a calculation of the kinetic 
The evaluation of the distorted ground-state energy contribution to the elastic energy only the one-particle 



density matrix 



i?2[u] is in principle straightforward if the exact one- 
particle and two-particle density matrices of the ground- 
state are known. In this section we focus on the con- 
struction of the kinetic contribution, which, as we will 
show, leads to a local equation of motion, which involves is needed. 

I 



p(r,r') = {^o|* + (r)*(r')|Vo) 



(49) 



The kinetic energy of the distorted state is 

T[u] = ±J dv^g^d^g-^g-y^p^v')]^, , 



(50) 



which reduces to the kinetic energy of the ground-state when u = 0. Expanding the above expression to second order 
in u we arrive, after some laborious algebra (see Appendix C for the derivation) to the following expression 



T 2 [u] = J dr^T^o 



+ -^(diu,u vl/ )(dnU 1/v ) 



where 



T, 



2m 



(d^ + dvd'J p(r,r') 



Am 



(51) 



(52) 



is the equilibrium stress tensor. Notice that T 2 [u] is a local functional of u, i.e. it presents no coupling between 
displacement fields at different positions. Taking the functional derivative with respect to u(r) we arrive at the 
desired expression for the kinetic force density: 



= d a [2T v „ u va +T va fidu,u v ] - ■^—d^dJnod^'V ■ u) 
ou^ 4m 

+ -^9 U {2(V 2 n )w^ + (9 i ,n )9 jU V • u + (d^n^d^X? ■ u - 2d^ [{d a n )u va ]} . 



(53) 



D. Calculation of 8E2[u]/8u — potential part 

To calculate the potential energy functional W[u] we need the two-particle density matrix of the ground-state: 

p 2 (r,r') = (V>o|* t (r)*V)%')*(r)lV'o> • (54) 
For a system of electrons interacting via Coulomb interaction (charge — e) we have 

P2(r,r') 



W[u] = %- f dr [ dr'- 

L J 2 J J |r + u 

Expanding to second order in u we easily obtain 



(r) - r' - u(r')| ' 
W 2 [u} = ~ Jdr [ dr'[u„(r)-u^r')]K^(r,r')[u„(r)-u v (r')} 



where 



K^(r,r') = p 2 (r,r')d fl d' v] 



r 



(55) 



(56) 



(57) 



Finally, taking the functional derivative with respect to 
u(r) we get 



Thus, the inclusion of interactions transforms our equa- 
tion of motion into an integro-differential equation. No- 



- = J dv'K^{v,v')[u v {v) u v {r> 



)]. (58) 
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tice, however, that the interaction contribution vanishes 
if u(r) is constant in space, as expected from the trans- 
lational invariance of the interaction. 



E. Equation of motion for one-dimensional systems 

The formulas presented in the preceding subsections 
simplify dramatically in one-dimensional systems, where 
the displacement field has only one component, u{x), the 
strain tensor reduces to the derivative of the displace- 
ment field u'(x), and the equilibrium kinetic stress tensor 
reduces to a scalar 



T (x) 



1 

77! 



d x d x >p(x,x')\ 



(59) 



Then the combination on the last line of Eq. vanishes 
and we are left with the simpler expression 

- = (37W)' - -±- (n u")" (60) 

OMp 4m 

where the primes denote derivatives with respect to x. 
The complete equation of motion for one dimensional 
systems is thus 

mnodtU — — uquVq + (3Tou)' — - — {n^u")" 

4m 

+ / dx K{x, x')[u(x) — u(x')] — noV{ , 



(61) 

where K(x, x') is given by the one-dimensional version of 
Eq. (f57|) . We will make use of this form of the equation 
of motion in the model applications presented below. 



III. CURRENT-DENSITY FUNCTIONAL 
APPROACH 

Our discussion thus far has not relied on time- 
dependent current density functional theory, except on a 
very abstract level, i.e. as a basis for the statement that 
the stress tensor must be a functional of the current den- 
sity. The formulas presented in the last two subsections 
relied on the knowledge of the exact density matrices 
p and p2 of the many-body ground-state - two quanti- 
ties that are amenable to treatment by powerful numer- 
ical techniques (e.g. the quantum Monte Carlo method) 
which have little in common with DFT. Before proceed- 
ing, we wish to clarify how the time-dependent CDFT 
can help us in more concrete ways when the exact p and 
pi are not known, which is by far the most common case. 

One of the main ideas of TDCDFT is that the current 
and density evolutions of the interacting many-body sys- 
tem can be simulated in a non-interacting many-body 
system subject to an effective time-dependent vector po- 
tential which includes the Hartree electrostatic potential 



and dynamical exchange-correlation (xc) effects. This 
non-interacting system is known as the "Kohn-Sham ref- 
erence system" and its ground-state density coincides 
with the exact ground-state density of the interacting 
system, i.e. no(r). The potential that produces this ex- 
act ground-state density in the Kohn-Sham reference sys- 
tem is known as the static Kohn-Sham potential and is 
usually written as 



K, (r) = Vb(r) + V H ,o(r) + V xc<0 (v) . 



(62) 



where Vh.q and V xc .q are, respectively, the Hartree poten- 
tial and the xc potential of the ground-state. These static 
potentials should not to be confused with the additional 
dynamical Hartree and xc potentials, which appear when 
the system is not in equilibrium. 

The idea is now to apply our continuum mechanics 
formulation directly to the Kohn-Sham reference sys- 
tem. There is a small technical problem in doing this, 
namely the time-dependent xc vector potential A. xc (r,t) 
that acts on the Kohn-Sham system has in general a 
transverse component, which cannot be represented as 
the gradient of a scalar potential. Indeed, a complete 
representation of A requires that we introduce both an 
electric field ~E XC and a magnetic field TS X c- The inclusion 
of the xc magnetic field does not create any difficulties 
in principle (see Ref. and leads to the appearance 
of a Lorentz-force term in the equation of motion for the 
current. But this Lorentz force term can be safely disre- 
garded in the linear response approximation, because it 
has the form j x B which is of second order in the devi- 
ation from equilibrium. Thus we can take into account 
dynamical xc effects simply by adding the force -eE IC 
to the driving force — noVVi on the right hand side of 
Eq. (fTS")) . All this considered, our equation of motion 
takes the form 



m<9 2 u(r,^) + (u • V)W S 



s,0 



1 ST s2 [u] 



n <5u(r) 

- V[Ti(r,t) + VH,i(r,t)]-eE ieCl i(r > t), (63) 

where Vr,i is the first-order term in the expansion of the 
time-dependent Hartree potential, and E^i is the first 
order term in the expansion of ~E XC in powers of u(r). 
The non- interacting kinetic force density — 5T s a[u]/(5u(r) 
is given by Eq. (|53l) in which, however, the equilibrium 
kinetic stress tensor is replaced by the correspond- 
ing quantity for the Kohn-Sham reference system, i.e 



2m 

-T-V 2 n (5^ 

4m 



(64) 



where tpe(r) are Kohn-Sham orbitals for the ground-state 
and the sum over £ runs over the occupied orbitals. 

Assuming that the Kohn-Sham ground-state has been 
obtained by one of the available approximations for V Sj o, 
the remaining problem is to find a suitable approximate 
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expression for E Ic .i. The "natural" approximation, in 
the present context, would be the high-frequency approx- 
imation, which expresses E lc .i as the functional deriva- 
tive of the exchange-correlation energy functional with 
respect to u. In practice, since the latter is not known, 
one has to rely on more or less uncontrolled approxima- 
tions, such as the high-frequency limit of the local den- 
sity approximation proposed in Refs i 19 ' 28 " — - see Eq. 
(117) of Ref. [H (a general discussion of these approxi- 
mation can be found in Ref. [3ll ) . This approximation is 
local both in space and time and is obtained by applying 
an instantaneous geometric deformation to the homoge- 
neous electron gas. In the second respect it is perfectly 
consistent with our elastic approximation for the nonin- 
teracting kinetic force, but we must keep in mind that 
the latter is fully nonlocal. 

Unfortunately, the local deformation approximation 
for the xc potential suffers, like all electron-gas based 
approximations, from a serious defect: it fails to cancel 
the unphysical self-interaction that is contained in the 
Hartree term. This makes it unsuitable for the treat- 
ment of strongly correlated system, where the spurious 
self-interaction energy can be very large. More accu- 
rate approximations^ 2 - do not suffer from this defect, but 
are more difficult to implement. Furthermore, such ap- 
proximations would do little (apart from fixing the self- 
interaction problem) to capture the physics of strongly 
correlated electrons. Alternatively, one could use func- 
tional explicitly designed for electronic systems in the 
strong coupling limit, such as the ones developed by 
Seidl, Perdew and Kurth^ and, more recently by Seidl, 
Gori-Giorgi, and Savin 34 . 



IV. THE CALCULATION OF EXCITATION 
ENERGIES 

A. The eigenvalue problem 

An immediate application of our equation of motion 
is the calculation of excitation energies^! 4 - To this end 
we turn off the external potential V\ and consider the 
homogeneous equation 



TOno(r)<9 t 2 u(r, t) = — 



5E 2 [u] 
<5u(r,i) 



(65) 



Fourier-transforming with respect to time and carrying 
out the indicated expansion of the energy to second order 
in u we get 



dr' 



S 2 E[u] 



5u^(r)8u v (r') 



u u lr ,w 



(66) 

Although this expression is not the most useful in prac- 
tice, it does bring forth some important features of the 
problem. First, because the kernel of the integral equa- 
tion is a symmetric second functional derivative, we are in 



the presence of an essentially hermitian eigenvalue prob- 
lem. More precisely, the problem is hermitian with re- 
spect to a scalar product with weight no(r), i.e. defined 
as 



(f,g) = / drn (r)f(r)g(r), 



(67) 



where / and g are arbitrary functions. This can be seen 
by rewriting Eq. (|66"T) as an equation for u = ^/no(r)u(r) 
and noting that this equation has the form of a standard 
eigenvalue problem with a symmetric kernel: 



dr'- 



1 



S 2 E[u] 



1 



,o y/mn (r') 



This means that all the eigenvalues will be real, and 
eigenfunctions u corresponding to different eigenvalues 
are orthogonal with respect to the ordinary scalar prod- 
uct. It follows that the original eigenfunctions u are or- 
thogonal with respect to the scalar product defined by 
Eq. (f6"T|). Second, the kernel of the integral equation is 
positive definite, because E[u] has an absolute minimum 
at u = 0, which corresponds to the ground-state energy. 
For this reason, it is guaranteed that all the eigenvalues 
are positive. The square roots of these eigenvalues are 
the approximate excitation energies of the system, start- 
ing from the ground-state. The eigenfunctions also have 
a simple interpretation as approximate matrix elements 
of the current density operator between the ground-state 
and the excited state under consideration, divided by the 
ground-state density. We will show this more clearly in 
the next section. 



B. Derivation from linear response theory 

Additional insight into the significance of the eigen- 
value problem for the excitation energies is obtained by 
deriving the equation of motion directly from the lin- 
ear response of the current density to an external vector 
potential in the high-frequency regime. To this end we 
write 

fr(r,w) = f dr' XliV (r,r', u)A Vil {v' ,w) , (69) 

where Xp„(r, r', w) is the current-current response func- 
tion, and notice that, at high frequency, this function has 
the well-known expansion^ 



X^(r,r',w) 



(70) 



where the first term (diamagnetic term) is frequency- 
independent and 

M„(r,r') = -m 2 (^o|[[^o,^(i-)],i(r')]l^o) , (71) 
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is the first spectral moment of the current-current re- 
sponse function. ([A-B] denotes the commutator of 
two operators A and B, and ipo is the undeformed 
ground-state wave function.) Now, replacing j(r, oj) = 



-iwno(r)u(r, ui) and Ai(r, w) 



VVi(r 



and solving 



for Wi to leading order in 1/w 2 we obtain 

d^Vi = - — / dr'M^v, r'K(r') . (72) 

n J 

This is equivalent to our equation of motion (14T>1) if and 
only if 



S 2 E[u] 



fe Al (r)(5u ly (r / ) 



(73) 



To show that this is indeed the case we observe that the 
deformed ground-state wave function can be expanded as 



(74) 



where ipQ is the undeformed ground-state wave function 

and V'o 2 ^ are corrections of first and second order 

in u respectively. The various corrections are not mu- 
tually independent. If ■f/'ol 11 ] is normalized to a constant 
independent of u, then we must have 



(M4 2) ) + i4 2) \i>o) = -(4 1) \4 1) ) 



(75) 



Taking this into account it is easy to verify that the sec- 
ond order correction to the energy is 



E 2 [u} = (4 1) \H -Ea\^ 



(76) 



where Eq is the ground-state energy. Finally, we observe 
that the first-order correction to the ground-state wave 
function is given by 



(i) 



im / drj(r) ■ u(r)\ip ) 



(77) 



where we have used the fact that the momentum density 
operator mj(r) is the generator of a local translation of 
all the particles in an infinitesimal volume located at r. 
Thus, the operator on the right hand side of Eq. ([77| per- 
forms different translations by vectors u(r) at different 
points in space, i.e. precisely deforms the ground-state 
according to the displacement field u(r). Substituting 

the above expression for \tpQ ) into Eq. (|76|) for ^[u] 
one can easily verify that 

E 2 [u] = lJ dr J dr'?v(r)AV(r,r'K(r') , (78) 

with M M1/ (r,r') given by Eq. ([71]). This establishes the 
validity of Eq. ([731). 



C. Eigenfunctions and sum rule 

The spectral representation of the kernel of our equa- 
tion of motion gives additional insight into the nature of 
our approximation and shows clearly where things can 
go wrong. Namely we can write 

M M „(r,r') = m 2 ^w„ {[j Al ] „(r)[j i ,] Il0 (r') 

n 

+ [>]o„(r')[^]«o(r)} , (79) 

where oj„o are the exact excitation energies of the system 
from the ground-state (0) to the n-th excited state, and 
lv]no( r ) = j M (r) 1 0) are the matrix elements of the cur- 
rent density operator between the corresponding states. 
An exact linear equation for the excitations would have 
to give ±oj„o as excitation energies and ^ 1 (r)[i/j]on(r) 
(= tt-q 1 (r) [j^] * (r) ) as the corresponding eigenfunctions. 

In general this will not be the case. However, for the 
special case of a one-particle system, in the absence of a 
magnetic field, we can show rather easily that [j M ]on(r) = 
— [j/i]no(r) and furthermore that n ( ^ 1 (r)[j At ]on(r) is indeed 
an eigenfunction of the operator n ( ^ 1 (r)M Miy (r, r') with 
eigenvalue muj 2 l0 . This follows from the orthonormality 
relation 



2m 



dr 



\jv]on(r)\j v ]ko(r) 
n (r) 



(80) 



which is valid for one-particle systems (in the absence of 
a magnetic field) and is proved in Appendix D. Then the 
eigenvalues of the equation of motion ([72"j) with V\ = 
are u = ±w„o as they should be. This result is per- 
fectly consistent with our previous observation that the 
"elastic approximation" is not an approximation at all 
when it comes to one-electron system, due to the lack of 
retardation effects in such systems. 

In general, in a many-particle system the matrix ele- 
ments of the current density operator between the ground 
state and different excited states are not necessarily or- 
thogonal. Indeed, we will see that two completely dif- 
ferent excited states can produce the same eigenfunction 
for the displacement field, up to a proportionality con- 
stant. The reason why this can happen is that the ex- 
act equation of motion for the displacement field of a 
many-body system is not an eigenvalue problem (even 
though it is linear) due to the frequency dependence of 
the kernel. As a result, the normalization of the solu- 
tions becomes relevant: two "eigenfunctions" that differ 
by a mere proportionality constant can result in different 
excitation energies when the kernel of the linear equa- 
tion is itself energy-dependent. In such cases, the elastic 
approximation will fail to resolve the different excitation 
energies, replacing them by a single excitation energy at 
an "average" value. 

In spite of this shortcoming, an exact sum rule can 
be established, which relates the exact eigenvalues u>\ 
of the elastic eigenvalue problem to the exact excitation 
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energies w n o- The sum rule reads: 



n 

where the "oscillator strengths" 



a, .2 

; A0 



/; 



_ 2m |/ rfrjo„(r) • u A (r)| 



(81) 



(82) 



and where u^(r) is the solution of the elastic eigen- 
value problem normalized with respect to the scalar 
product ([67|) . ui\ is its eigenvalue, and the sum runs 
over all the exact eigenfunctions. Further, the oscillator 
strengths satisfy the sum rule 



(83) 



for all A. The proof of these results is presented in Ap- 
pendix E. 

From this vantage point we see that the elastic ap- 
proximation is the extension of the well-known collective 
approximation^ of the homogeneous electron gas to in- 
homogeneous systems. Each eigenvalue uj\ of the elastic 
equation of motion is a weighted average of exact exci- 
tation energies, with a weight controlled by the overlap 
of the exact current matrix element with the eigenmode 
u^(r). Since u,x(r) form an orthogonal basis in the space 
of displacements, we can say that uj\ represents the av- 
erage energy of exact excitations in the "direction" A. 
Thus, the full excitation spectrum is replaced by a set of 
spectral lines, one for each orthogonal direction in dis- 
placement space, and each one carrying the entire and 
exact spectral weight for that particular direction. 



characterized by a wave vector q. Of these there are 
two kinds: longitudinal, in which u is parallel to q, and 
transverse, in which u is perpendicular to q. The expres- 
sion (fS"3")) for the kinetic force density reduces to 



<5T 2 [ul 2 .... , . 9 , no 2 . 

— iU = _ nt (n 2q q-u )+q 2 u + -*-q(q. u . 
ou 6 4m 



(86) 



The force density from potential energy, Eq. (|5"5|) . is given 
by 



SW 2 [u} 



[K^(0) - K^c^uv , (87) 



where 



where v(q) = 4ire 2 /q 2 is the Fourier transform of the 
Coulomb potential in three dimensions, and /?2(q) ~~ the 
Fourier transform of p 2 (r — r') - is related to the static 
structure factor S(q) in the following manner^ 



p 2 (q) = n[S(q) - 1] . 



(89) 



Thus, we get 

_ 6W 2 [u] = 
Su 



-n J dq'[5(q - q') - S(q[)W)q[[q! ■ u(q)] 



V. MODEL APPLICATIONS 

For orientation we now examine the application of our 
theory to a few simple models. 

A. Homogeneous electron gas 

In a homogeneous electron gas the ground-state den- 
sity no — n is independent of position. The equilibrium 
kinetic stress tensor has a constant value 

2 

Tn„,o = 2 n *( n )<W > ( 84 ) 

where t(n) is the kinetic energy per particle. The two 
particle density matrix ^(r, r') is a function of |r — r'|. 
In such a homogeneous system the displacement eigen- 
functions are simply plane waves 



Finally, in order to take into account the neutralizing 
background of positive charge (required for the stability 
of the electron gas), we add the external potential 

y (r) = ^(r-q) 2 , (91) 

where w 2 = 47rne 2 /m is the square of the plasmon fre- 
quency and q is the unit vector in the direction of q. 
Notice that this potential is assumed to vary only in the 
direction of q, because it is only in this direction that the 
displacement generates boundary charges: the system re- 
mains perfectly homogeneous in the direction perpendic- 
ular to q4^ The corresponding force in the equation of 
motion (fTSj) is 

u- Va,y = m^(q.u)g (1 . (92) 



u(r,w) = u(q,w)e 4(q - r -" t) 



(85) 
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Putting everything together, the equation of motion takes the form 

m^u = -t(n) [2q(q • u) + q 2 u] + |^q(q • u) + m^(q . fi)q + — * J [S(q - q') - S(q')]q'[q' • u(q)] , 

(93) 

This can be further decoupled into longitudinal and transverse components denoted by ul and ut respectively. The 
corresponding eigenvalues are 

= ^ + 2t(n)^ + ^ + ^ f-p-iH- q'f [5(q - q') - 5(q')] , (94) 
^ to 4m z n J (ZTry 

and 

4(9) = ^(n)^ + ^ | ^_ (q x q f [5(q _ q /) _ 5(q /)] . (95 ) 



Of course, this is exactly what one would have obtained 
by assuming that the spectrum of the current-current 
response function consist of a single <5-function peak at 
lul or lot and requiring satisfaction of the first moment 
sum rule (see Ref. HU, Eq. 3.191 and Ref. [30h . 



In Fig. (fTJ we plot the excitation spectrum of the ho- 
mogeneous electron gas calculated from Eqs. (|94l) and 
([95]) . We have used the static structure factor calculated 
in Ref. by the quantum Monte Carlo method, and the 
kinetic energy has been computed from the parametrized 
correlation energy of Appendix B of Ref. l36l using the 
virial theorem. In the longitudinal channel, the exact 
spectrum is dominated at small q by the plasmon and 
at large q by free particle excitations (energy q 2 /2m). 
There are also also electron-hole pair excitations at lower 
q and ui, as well as multiple electron- hole pair excitations 
all over the plane. The elastic approximation replaces 
this complex spectrum by a single branch of longitudinal 
excitations which has the correct spectral moment. In 
particular, we get the correct dispersion of the plasmon 
at small q and the correct free particle behavior at large 
q. In the transverse channel the plasmon and the high-q 
free particle excitations are absent. The exact transverse 
spectrum consists primarily of low-energy electron-hole 
pair excitations. The current vector, k + q/2, where k 
and k + q are the momenta of the hole and the elec- 
tron, respectively, is essentially perpendicular to q when 
k and k + q both lie near the Fermi surface. On the 
other hand, high-energy excitations, with energy q 2 /2m 
are essentially longitudinal because the current vector is 
essentially parallel to q when q is much larger than the 
Fermi momentum. This is consistent with the fact that 
the frequency of our transverse collective mode u>t(<?) 
grows linearly with q at large q. 



B. Linear harmonic oscillator and hydrogen atom 

In order to demonstrate the exactness of the formu- 
lation for one-electron systems let us now consider the 
canonical examples of the one-dimensional harmonic os- 
cillator and the hydrogen atom. 

For a harmonic oscillator of natural frequency luq, ex- 
ternal potential Vq(x) = muj^x 2 /2, and equilibrium den- 
sity Uo(x) = 6 ^ £ , where I = (mwo) -1 ^ 2 , the equation 
of motion ([STj) reduces to 

\u"" - xu'" + (x 2 - 2)u" + 3xu' - ^—^u = . (96) 
4 u 2 

Solving the eigenvalue problem with the boundary con- 

1/2 

dition of n Q (x)u(x) — > as |x| — > oo, we obtain the 
exact excitation spectra oj n — ±nQ, where n — 1,2,... 
The corresponding cigenfunctions are 

Wn(x) oc H n -i(x) (97) 

which are mutually orthogonal with respect to the scalar 
product (|67[) . These are indeed proportional to the ma- 
trix elements of the current density operator between the 
ground-state and the n-th excited state. 

A similar calculation can be done for hydrogen-like 
atoms of atomic number Z. Focusing for simplicity on 
excitations of spherical symmetry, we introduce a radial 
displacement field u r (r) which depends only on the radial 
coordinate r. Then u r satisfies the equation 

+ ^<-(^ + !l)^ = ' (98) 

where the primes now denote derivatives with respect to 
r. Solving the eigenvalue problem with boundary con- 

1/2 

dition n ' (r)u r {r) — > for r — > oo yields the correct 
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FIG. 1: (Color online) Longitudinal (a) and transverse (b) 
modes for a homogeneous electron gas at r s — 1,3,5. Wave 
vector q is in units of fcj and frequency is in units of 2Ef- 
The curves labelled uj+ (q) and ui- (q) are the boundaries of the 
electron-hole continuum.—. Single particle excitations exist 
for cj—(q) < u < Lj+(q), whereas multiparticle excitations 
are distributed all over the plane. The elastic approximation 
replaces the exact spectrum by the two branches uih{q) and 
u)t(<i), which carry the entire spectral weight. Notice that 
the r 3 -dependence is barely discernible on a large q scale, but 
becomes clearly visible at smaller q as shown in the insets. 



excitation energies oj n — (Z 2 /2)(l — 1/n 2 ) (n = 1,2,..). 
The corresponding eigenfunctions are given by Laguerre 
polynomials: 



u n (r) oc L 2 n _ 



(99) 



relative variable, and analytically in the limit of strong 
correlation. The hamiltonian is 



Ho = — 

4m 



2 XT' 2 Jr 2 2 

mcUnX H 1 cu n x 

' m 4 



VP 



:, (100) 



where X = Xl + X2 and P = Pi + P2 are, respectively, 
the coordinate and the momentum of the center of mass 
and x = xi — X2, p = Pl ~ P2 are the coordinate and the 
momentum in the relative channel. Notice that for a 
fixed strength e 2 of the interaction we can go from the 
weakly correlated to the strongly correlated regime by 
varying the value of ujq: ujq — > oo corresponds to the non- 
interacting limit, and uo — » to the strongly correlated 
limit. 

The ground-state wave function is 



^o(x 1 ,x 2 ) = 1p (X)(j) (x) 



(101) 



where ipo(X) and 4>o(x) are, respectively, the ground- 
state wave functions of the center of mass and of the 
relative hamiltonian. The general excited state is 



^ nm (x 1 ,x 2 ) 



(102) 



where tp n (X) is the wave function of the n-th excited 
state of the center of mass hamiltonian and <j) m {x) is the 
wave function of the m-th excited state of the relative 
hamiltonian. 

The ground-state of the system is a spin singlet (S = 
0) and for this reason in the following we consider only 
singlet states, which are connected to the ground-state by 
the current density operator. The relative wave function 
for such states is symmetric: (f> m (— x) — (f> m (x). This 
wave function has 2m nodes, of which m with x > and 
m (symmetrically placed) with x < 0. The center of mass 
wave function, 



MX) = H n 



X 



-X 2 /2 



(103) 



can be either symmetric or antisymmetric, depending on 
the parity of n, and has n nodes. Here £ cm = \Jh/2rauiQ. 
The ground-state has n = 0, m = and all the other 
states are characterized by positive values of the integers 
n and m. 

In the non interacting limit (ojq — > oo, or e 2 — > 0) the 
relative wave function is 



C. Two-electron systems 

As a final example let us consider the case of two 
electrons repelling each other with the "soft" Coulomb 

2 

potential . e in a one-dimensional parabolic 

trap of natural frequency ljq (the cutoff a > serves 
to eliminate the pathological behavior of the interaction 
at x 1 = x 2 ). This is a model that can be solved nu- 
merically thanks to the separation of center of mass and 



4>m(x) = H 2m 



-x 2 /2ll 



(104) 



where £o — ^/2fijmu)Q. The ground-state density is a 
gaussian centered at the origin. The excitation energies, 
expressed in units of uiq, are sums of the excitation ener- 
gies of two identical harmonic oscillators: 



lim 

u — >oo ujq 



2m 



(105) 
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FIG. 2: (Color online) Evolution of the excitation energies 
for two electrons in a one-dimensional harmonic trap. Solid 
lines denote exact excitation energies, labelled by (n, m) as 
explained in the text. Solid dots indicate the calculated eigen- 
values of the QCM equation of motion. Crosses on the right 
denote the strong-correlation limit of the eigenvalues, given 
by Eq. UTT3J. 



with n, m non- negative integers. The degeneracy of the 
excited states is the number of integers less or equal n + 
2m with the same parity as n + 2m, i.e. 



2m 



(106) 



where [y] denotes the integer part of y. The displacement 
field associated with the (n, m) excitation is 



i(x) oc H n+2m -i{x/\o) , 



(107) 



where Ao = \J l\ m + {Iq/2) 2 — y/h/muJo. The parity is 
j_ljn+2m-i anc j ^ e num b er f nodes in n + 2m — 1. 

The situation is quite different in the strongly corre- 
lated limit (wo — > 0, or e 2 — > oo). The relative Hamilto- 
nian reduces to a harmonic oscillator of frequency ojqVS 
with equilibrium distance so = (2e 2 /majg) 1 / 3 (we assume 
a <C xq). The ground-state wave function for the rel- 
ative motion is a symmetric linear combination of two 
Gaussians of width — (2S/'\/3mwo) 1 ' 2 centered at 
x\ — X2 = ±xq. The corresponding ground state den- 
sity tiq{x) consists of two Gaussian peaks of the width 

Aoo = ^JP cm + (4o/2) 2 = [h(l + ^^VSmiOo] 1 / 2 cen- 
tered at x = ±xo/2. The excitation spectrum has the 
form 



and the degeneracy is completely removed. The displace- 
ment field is analytically found to be 

(x) cx if n+m _i \ X ° ) 9(x) 



lim = n 



iv3, 



(108) 



(-11'"//,.-,,, , ( X \ X0lT )ti(-r)m\)) 



The parity is (—1)" 1 and the number of nodes is 2(n + 
m— 1) + mod(n — 1,2), where mod(n— 1.2) = n — 1 (mod 
2). 

The evolution of the lowest-lying energy levels with 
given value of the pair n, m as a function of wo is shown 
by the solid lines in Fig. @. Some of the displacement 
fields of the low-lying excitations in the strongly corre- 
lated regime (Eq. (|109l) ) are shown in Fig. 

From these figures we see that the displacement field of 
the (1,0) excitation, which corresponds to a rigid trans- 
lation of the center of mass, is uniform in space, while 
the displacement field of the (0, 1) excitation, which cor- 
responds to the classical breathing mode, changes sign 
around the origin. The (1, 0) and (0, 1) excitations corre- 
spond to the classical phonon modes of a system of two 
localized particles. The remaining excitations are quan- 
tum mechanical in character, as can be surmised from the 
fact that their displacement fields (in the strongly corre- 
lated regime) have significant variation over the regions 
where the density has peaks, i.e. the places where the 
particles would be classically localized. These modes de- 
scribe the dynamics of the wave function of the localized 
electrons. 

Looking at the figures we observe that, in the 
strongly correlated regime, there are groups of excited 
states (e.g. {(0, 2), (2, 0)}; {(3, 0), (1, 2)}; {(2, 1), (0, 3)}; 
{(4,0), (0,2), (0,4)}, such that all the excited states 
within one group produce the same displacement field, 
up to a normalization constant. In general, states with a 
given value of n + m and the same parity of m have the 
same displacement field, but different energies. Clearly, 
this is a feature of the exact solution that cannot be 
reproduced by any linear eigenvalue problem with a 
frequency-independent kernel. 

The phenomenon of different excited states produc- 
ing the same displacement field occurs also in the non- 
interacting limit: all the states with the same value of 
n + 2m (e.g. (0,1) and (2,0)) have the same displace- 
ment field. But, in this case, the states with the same 
displacement field also have the same energy: therefore 
the noninteracting excitation energies can be accurately 
reproduced by a linear eigenvalue problem. 

Let us now see what our elastic equation of motion (|6ip 
predicts for this system. The kinetic part of the equilib- 
rium stress tensor Tq(x) works out to be 
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(0,1) o 
(0.2), (2,0) • 
(2.1), (0,3) - 



'OOOOOOOOOOO 




FIG. 3: (Color online) Top panel: the displacement field 
u„ m (x) for (n, m) = (0, 1), (0, 2) and (0, 3) in the strong corre- 
lation limit. Bottom panel: the same for (n, m) = (1, 1), (1, 2) 
and (2, 1). The thin solid lines represent the density profile. 
The large value of the displacement field for x ~ does not 
have a physical significance since the density is exponentially 
small in that region. 



T (x) = — dy 
m 



x + y 



L l A 1 \ /' ( X + V 



x + y 



(110) 



where ipo an d 4>o are the ground-state wave function in the center of mass and relative channel respectively. The 
interaction kernel K(x, x') is given by 



K(x,x') 



-2{x-x'Y+a 2 
[(x - x 1 ) 2 + a 2 ] 5 / 2 



(p (x - x')ip 



x + x' 



(III) 



We now have all the input that is necessary to set up and solve the fourth-order integro-differential equation (|61[) . 
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In the limit of weak correlation (cjo — > oo) the eigen- 
values of the intcgro-diffcrcntial equation coincide with 
the exact (degenerate) excitation energies. This is un- 
derstandable, since in this limit the two electrons are de- 
coupled and the excitation spectrum of the two-electron 
system coincides with that of a single electron starting 
from its own ground-state. This spectrum, as we have 
seen, is exactly reproduced by our equation of motion. 
Unfortunately, this nice feature of the present model can- 
not be extrapolated to general systems. If, for example, 
the system contains more than two electrons, then even 
in the non-interacting limit a generic excitation will en- 
tail the transition of a single electron from an occupied 
orbital that is not the ground-state orbital to an unoccu- 
pied one. Such an excitation will not be described exactly 
by our method. 

In the limit of strong correlation our integro- 
diffcrcntial equation can be solved analytically, as shown 
in Appendix F. The eigenfunctions can be classified as 
even (+) or odd (-) and are given by 



Even modes 



Uk,±{x) oc H k ( I 9(x) 



-x), (112) 



where k is a non- negative integer. The number of nodes is 
2k for even eigenfunctions, 2fe + l for odd eigenfunctions. 
The corresponding eigenvalues are given by 



lim 

uo— >o ujq 



2 + 3V3fc + 6k(k - 1)(2 - VS) 
T (-l) k (2- V3) k } 1/2 . (113) 



Notice that, within each symmetry sector (even or odd), 
the eigenvalues increase monotonically with increasing k. 

In Table I, fourth and fifth column, we present a de- 
tailed comparison between the exact excitation energies 
and the eigenvalues of our equation of motion in the 
strong correlation limit. For the sake of clarity, we list 
the excitations that produce even displacement fields and 
those that produce odd displacement fields separately. 
The elastic equation of motion can also be solved numer- 
ically, and the results are in very good agreement with 
the analytical solution. This, and the fact that the sum 
rule ([51]) is satisfied with good accuracy, builds our con- 
fidence in the numerical solution. 

In Fig. ^ we present the numerical results for some of 
the lowest-lying excitations as a function of ujq. We can 
immediately see that the "non-degenerate" excitations, 
by which we mean the excitations (1, 0) (0, 1), and (1, 1), 
which arc uniquely associated to a given displacement 
field, are rather well reproduced by our calculation for 
all values of ojq . On the other hand, the "degenerate 
excitations" , which yield the same displacement field but 
have different energies, are replaced by a single excitation 
of an average energy, in such a way that the total spec- 
tral strength of the group is preserved. Two examples of 
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TABLE I: Exact excitation energies E nm of the two-electron 
model with hamiltonian specified in Eq. (| lOOf) in the non- 
interacting limit (E^ m , second column) and in the strongly 
correlated limit (E^ m , fourth column). The eigenvalues of 
the QCM equations of motion (.E£+m-i,±) m the strongly 
correlated regime are listed in the fifth column. The top half 
of the table lists excitations with even displacement fields and 
the bottom half lists excitations with odd displacement fields. 
The third and the last column on the right list the number 
of nodes in the displacement field in the non interacting limit 
(iV ) and in the strongly correlated limit (N°°). 



this phenomenon are evident in Fig. ([2]): the (2,0) and 
(0, 2) excitations, which in the strong correlation limit 
(ujq —} 0) have energies 2luq and 3.464a>o respectively, 
are replaced by a single excitation - the fourth one in 
Fig. |J5) - which tends to the "average" energy 2.632wo- 
Similarly, the (3,0) and (1,2) excitations, which, in the 
ujq — > limit tend to 3wo and 4.464u;o respectively are 
replaced by a single excitation - the fifth one in Fig. @ 
- which tends to the "average" energy 3.942cjo- The pat- 
tern recurs for more complex multiplets of excitations, 
involving three or more states with the same displace- 
ment field and different energies. 

We notice that the displacement field associated with, 
say, the (n, m) excited state has a number of nodes that 
generally grows from n + 2m — 1 in the weak coupling 
limit to 2(n + m— 1) (odd n) or 2(n + m— 1) + 1 (even n) 
in the strong coupling limit. This effect is particularly 
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FIG. 4: (Color online) Evolution of the displacement fields 
of excitations (1,1), (3,0) (even) and (2,0), (2,1) (odd) as a 

— 1/2 

function of correlation strength \ = > as shown in the 

top left panel. Notice the variation in the number of nodes as 
X increases from the weakly correlated to the strongly corre- 
lated limit. 



pronounced for states of small m, and is absent in the 
n = states. Fig. ([JJ shows the evolution of the dis- 
placement field for the even excitations (1, 1) and (3,0) 
and for the odd excitations (2, 0) and (2, 1). We see that 
in the "non-degenerate" (1, 1) state, the number of nodes 
stays constant and equal to 2 as one goes from the weakly 
correlated to the strongly correlated regime. In the (3, 0) 
state the number of nodes grows from 2 to 4 nodes, so 
that the displacement field of this state becomes pro- 
portional, in the strong-correlation limit, to that of the 
much higher in energy (0, 3) state. The same behavior is 
observed in state (2,0), for which the number of nodes 
grows from 1 to 3, and in state (2, 1), for which it grows 
from 3 to 5. By this mechanism, states of very different 
energy end up sharing the same displacement field (up 
to a proportionality constant) in the strong correlation 
limit. 

Our discussion has been limited to singlet states (sym- 
metric wave function in the relative channel). It would 
be easy to extend the calculation to include triplet states. 
To this end, we simply replace the density, kinetic energy 
density, and pair correlation function of the ground-state 
(a singlet) by the same quantities calculated from the 
ground state in the triplet (S — 1) sector of the Hilbert 
space. The relative wave functions of these states are an- 
tisymmetric. The correct symmetry of the wave function 
is automatically taken into account through the ground- 
state properties, and does not further appear in the elas- 
tic equation of motion. 



VI. DISCUSSION AND SUMMARY 



The elastic approximation is, in a very precise 
sense, the extension of the well-known collective 
approximatio n 39 ' 40 of the homogeneous electron gas to 
non-homogeneous electronic systems. In the case of two 
electrons interacting by Coulomb potential in a harmonic 
trap, we have seen that the elastic approximation re- 
places groups of excitations characterized by the same 
displacement field by a single excitation that carries the 
oscillator strength of the whole group. In more com- 
plex systems, we do not expect to be able to identify 
small groups of excitations that share the same displace- 
ment field. All that can be said is that the displace- 
ments associated with different excitations will not be 
linearly independent. Each eigenfunction of the elastic 
equation of motion will overlap with many different ex- 
citations. However, the integrated spectral strength of 
the elastic eigenmodes will still add up to the correct 
value. For this reason, our approximation should be use- 
ful in dealing with collective effects which depend on the 
integrated strength of the excitation spectrum, such as 
the dipolar fluctuations that are responsible for van-der 
Waals attraction ! 37 ! 38 Other possible applications include 
possible nonlocal refinements of the plasmon pole ap- 
proximation in GW theory 41 - and studying the dynam- 
ics of strongly correlated systems, which are dominated 
by a collective response. As a byproduct we got an ex- 
plicit analytic representation of the exact xc kernel in 
the high-frequency (anti-adiabatic) limit42 This kernel 
should help us to study an importance of the space and 
time nonlocalities in the KS formulation of TD(C)DFT. 
It would be particularly interesting to try and interpo- 
late between the adiabatic and anti-adiabatic extremes to 
construct a reasonable frequency-dependent functional. 

The elastic equation of motion derived in this paper 
relied on the knowledge of the exact density matrices 
and of the ground-state. In many cases, these 
ground-state properties can be extracted from Quantum 
Monte Carlo calculations. When this cannot be done, one 
can still resort to density functional theory, i.e. apply the 
QCM formulation directly to the Kohn-Sham system, in 
which case we do not need the exact ground-state den- 
sity matrices, but only the ground-state KS orbitals and 
a reasonable approximation for the exchange-correlation 
field. While the standard KS method treats the nonin- 
teracting kinetic stress tensor exactly, our method should 
be computationally more agile, for large systems, since 
it does not involve time-dependent orbitals and/or the 
inversion of large linear response matrices. 

It remains a challenge to extend the present formal- 
ism to the nonlinear regime, as well as including external 
magnetic fields and spin-orbit interactions. 
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Appendix A: Stress tensor operator 

From the evaluation of Eq. ([21)1) at the Euclidean met- 
rics 5 M „ = Sfnv we gefci&±2 



(Al) 



where 



2m 



(A2) 



and 



x [ dA/5 2 (r + Ar',r- (1 - A)r'). (A3) 
Jo 

Here ^(r) is the field operator, 

p 2 [r, r') = * f (r)* V)*( r ')*( r ) (A4) 

is the diagonal two-particle density operator, and w(r) is 
the interaction potential. 



Note that the identity relates the functional derivative 
of S with respect to the displacement to the functional 
derivative with respect to the metric/deformation tensor, 
which physically means a connection of the force to the 
stress. 

To prove Eq. (|B2p we consider a small variation of the 
function r^(£): ?>(£) h- » r^(^) +Sr fl (£). The correspond- 
ing variation of the functional S takes the form 

SS= [ dt^-Sr M (e) (B3) 
J or ^ 

On other hand, the variation of r M (£) induces the follow- 
ing variation of the metric tensor: g^(^) > gfu>(£) + 
$9nu(£), where 



d8r a dr a dr a dSr a 



Hence the variation of S can be also written as 



(B4) 



SS = [ de-^Sg^it) 



,^ SS f d5r a dr a | dr a dSr a . 



Performing the partial integration in the right hand side 
of Eq. (jB5[) . and using the symmetry of the tensor g^ 
we reduce Eq. (jB5[) to the following form 



SS 



M£) ( B6 ) 



ss = - 



The direct comparison of Eqs. (|B3[) and (jB6[) proves the 
announced identity (|B2|) . 

Finally, to make a connection to Eqs. (13"6l and (f3"5| we 
set S[g[j_v] — H u [9im>] an d take the expectation value of 
Eq. (|B2|) in the state \ip(t)). Then, linearizing with re- 
spect the displacement u(£), we find that the right hand 
side of Eq. ()B2|) becomes identical to the right hand side 
of Eq. <J36J), while the left hand side of Eq. J§2j is exactly 
equal to the right hand side of Eq. 



Appendix B: Derivation of the force identity 
Eq. flggj) 

In this appendix we derive an identity which is used 
in Sec. II A to identify the right hand sides of Eqs. (|56")l 
and ([3"8"1) . Namely, we consider a functional S[g^] of the 
following metric tensor 

(IV (IT 

9^t) = Zr^Ti r a ($)=£ a + u a ($) (Bl) 
and prove that the following equality holds 

or ,t d£ a \d£/3 dg aP J d£, a \dfy J 

(B2) 



Appendix C: Derivation of Eq.(|5ip 

In this appendix we derive the linearized form of the 
kinetic energy Eqs. ([ST]) for the instantaneously distorted 
ground state. 

We start with the general nonlinear expression for 
the kinetic energy T[u] in the elastic approximation [see 
Eq. (USD] 

T[U] = i / d ^9» v dM- l l\r)g-V\Op{£A')]t=e, 

( C1 ) 

where p(£, £') is the exact ground state one-particle den- 
sity matrix, and g^iC) and g(£) are, respectively, the 
inverse and the determinant of the metric tensor g^i^) 
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that is the functional of the displacement u(£), which is 
defined by Eq. (|BTt 

Our aim is to expand the functional of Eq. (|Clj) to 
the second order in the displacement field, i. e., to the 
first non- vanishing contribution corresponding to the lin- 
earized theory. 

First we explicitly calculate the derivatives in the right 
hand side of Eq. (JClJ and set £' = £. As a result Eq. (1UT) 
reduces to the form 

T[u] = J cLW{K^ + |?-(cklnys)(cU n V?) 

- ^[(^lnVg)(a^ ) + (^lnV5)(^n )]}, (C2) 

where n.o(£) — is the ground state density, and 

M0 = ^[WOl«-. (C3) 

Making use of the following representation for ^/g, 

^ = dBt (fe)' 
we can write the derivative of In y/g as follows: 

a„ tav5 = a „ lnd et(^)=§^|r <C4) 



It is now straightforward to expand the right hand side 
of Eq. (|C4[) to the second order in u: 



In yfg « d^daUa - {dpu^d^daUp . (C5) 

Next we consider the covariant tensor g^ v (the inverse of 
9y.u)- 

9^ = [Sp, + d^u v + d v Uft + (9 M M a )(9„« a )] 1 . (C6) 



Expanding the inverse matrix in Eq. (|C6[) to the second 
order and expressing the result in terms of the strain 
tensor u M „ we get 



g^ v ~ 5^ - 2u t _ iu + Au m u va - (d fl u a )(d u u a ) . (C7) 



Finally, we substitute Eqs. (|C5j) and jC7]) into Eq. (|C2|) 
and keep terms up to the second order in u. The second 
order contribution to T takes the form 



+ -^(^Ua)^^^)} . (C8) 

The last term in Eq. (|C8[) can be identically represented as follows 

^r^~(d^u a )d„(d a u^) = - \ n ° [^u^u^ - (d f _ l u a )(d fl u a )} - ^^-u^d^u^ . (C9) 
4m 8m 2m 

Note that the coefficient in front of the first term in Eq. (|C9|) and the corresponding coefficient in first term in Eq. (|C8[) 
are naturally combined into the kinetic stress tensor T^ q — 2K fJ/U — (5 AI1/ V 2 n /4m. Hence, inserting Eq. (|C9I) into 
Eq. (|C8p and integrating by parts terms proportional to <9„7io/2m we arrive at the following final representation for 
the linearized kinetic energy of the distorted state: 



T 9 



di^T^^Au^Uva - {d^u a )(d u u a )\ + ^^{d ll u aa ){d^u ul/ ) + ^[d^Uyad^U^a - d^Uv^d^Uaa] , (CIO) 



which is identical to Eq. ([5T]) . It is worth noting a con- 
venient feature of this representation - the last term in 
Eq. (IC10|) vanishes in all ID systems and for homoge- 
neous systems with tiq — const in any number of dimen- 
sions. 



Appendix D: Proof of the orthonormality relation 
for single-particle transition currents 

In this appendix we prove Eq. (|80| for one-particle 
systems in the absence of a magnetic field. 

The matrix elements of the current density operator is 

j „(r) = (0|j(r)|n) = ~ ^[V-oWn - , (Dl) 

where ipo,..., ip n are orthonormal eigenfunctions of the 
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one-electron hamiltonian, which can be assumed to be 
real if there is no magnetic field. Now observe that 



J0n(r) I V'oVl/'n - 



2m 



. — <0 O V ( — 
2m V ip 



(D2) 

and jno(r) = — jon(r). Also, from the continuity equation 
we get 

V • (ipo'Vip n - ip n Vipo) = -2muj n Qip ip n . (D3) 
Combining these two equations we get 



-2mu n oipai>n ■ (D4) 



From this we see that 
dr 



2m jon(r) \ / 1 2m j fc0 (r) 



1 1 



2m y/uinoUko 
1 1 



dr 



2m y/uJnQLOko J ip 



to 



WfcO 



fink 



dri/j n ijj k 



(D5) 



which proves Eq. (|80l) . 



Appendix E: Proof of the sum rules (fglj) - ([83]) 

Our starting point is the first moment sum rule for 
the current-current response function (or "third-moment 
sum rule" for the density-density response function), 
which states that 

1 f°° 

-- / dww9mx^(r, r',w) = V" w„o[j AI (r)]on[>(r')]nO • 

« 

(El) 



On the other hand, the equation of motion (|72|) for 
the current density in the elastic approximation can be 
rewritten as 



dr 



' [u>H^S{v r') - AV(r, r')} J^lj^) 

(E2) 



where = d^Vx/tyuS), j v — —iumoU v , and 
M^(r,r') = p=L== M M1/ (r,r')- ' 



v /mn (r') 



X) w ™Q{bV( r )]onbV(r')] 



no(r) 



/o 



is a hermitian positive definite operator, which admits 
a complete set of orthonormal eigenfunctions. Let us 
denote by u A (r) these eigenfunctions, and by u\ their 
eigenvalues. The orthonormality relation reads 

J dru A (r) ■ UA'(r) = 8\\> 

and the completeness relation is 

^[u A (r)] M [u A (r')], = <5(r-r')V. 



(E4) 



(E5) 



The kernel itself can be written as 

M^(r,r') = ^^[u A (r)] M [u A (r')],. 



(E6) 



The equation of motion for the current density can be 
rewritten as 

J dr'( W 2 -^)[ UA (r)],[u A (r')],^^J,(r') 



W^Mr), 



(E7) 



Its solution is obtained by projecting both sides of the 
equation along the eigenvector u A . We get 



m 

J 

n 



(E8) 



where the subscript A denotes projection along u A . From 
this we obtain 



n {r) 



W 2 — V TO 



w a ,/ n o(r) r ~ 



[fiV)W^.(r'). 



(E9) 



Hence, the current-current response function in the elas- 
tic approximation is 

X ^(r,r', W ) = ^^(r-r') 



or — w 



^ ■/ n ° (r) ^(r)] Al [u A (O], 1 / n0(r,) 



Hi 



Hi 



(E10) 



Evaluating the sum rule we obtain 



1 



ID 



™o(r) r ~ 



fl A (r)] M [u A (r / )] l 



no(r') 



/ 777, i 

+ ^(r')]o„fc(r)]„o}^^y- (E3) -^(r)AV(r,r')^(rO 



(Ell) 
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On the other hand, the first moment sum rule (|E1[) can 
be rewritten as 



1 f°° 
n Jo 



^ % /^)M^(r ) r')Vn (r') . (E12) 

Comparing the last two equations we conclude that 

1 f°° 
7T Jo 

— / duw9fmx«(ry.w), (E13) 

i.e. the sum rule is satisfied in the elastic approximation. 

More pointedly, making use of Eq. (IE3|) the sum rule 
can be written in the form 



^w^[u A (r)] p [uA(r')]^ = 2m^o; n o 



from which it follows that 



— E w «o/n 1 



[j M ( r )]o" [>( r ')]»0 
\Aio( r ) i/n (r') ' 
(E14) 



(E15) 



where the "oscillator strengths" /„ are positive quantities 
defined as 



f A = 

J n 



2m\F, 



A 12 



W n 



with 



= y dru A (r) 



[j( r )]o» 
\A*o(r) 



As a final step we prove that 



(E16) 



(E17) 



(E18) 



for all A. This is of course nothing but the /-sum rule 

duj £ = — d^d(r-r ) , (E19) 

w Am 



which is manifestly satisfied by xiu virtue of the com- 
pleteness relation (|E5[) for Ua(i"). When applied to the 
exact response function the /-sum rule implies that 



Then we see that 



E 



^(rjW^vCrQU [uA(r)] M [uaQ')^ 
0J n y/no(r) \AoIr 7 ) 



dr / dr'n (r)S^6(r - r') 



/JuaWI,, [uA(r')]^ 



V«o(r) \/n (r') 

(E21) 



QED. 

Appendix F: Analytic solution of the elastic 
eigenvalue problem in the strong-correlation regime 



In this appendix we present an asymptotically exact 
solution of the ID continuum mechanics eigenvalue prob- 
lem for two particles confined by the harmonic potential 
Vo = ^itiuqX 2 and interacting with a soft-Coulomb po- 
tential 



w(x — x) = 



^/{x - x') 2 + a 2 



(Fl) 



At the exact many-body theory level the system is de- 
scribed by the Hamiltonian of Eq. (|100[) . Within our con- 
tinuum mechanics the excitation energies are obtained 
from the solution of the following "elastic" eigenvalue 
problem for the displacement u(x,t) — u(x)e~ lU!t (see 
Eq. (pjTj) in the main text) 



muj 2 n u(x) = —d 2 x [n Q d 2 x u{x)\ - 3d x [T d x u(x)] + mu) 2 n u(x) +2 1 dx'[d 2 x w{x - a/)]#o(x, x')[u{x) - u(x')] . (F2) 



Am 

Here ^o(x,x') is the ground state two-particle wave function, which in this case coincides with the square root of 
the two-particle density matrix, no(x) is the ground state density, and T (x) is the ground state kinetic stress tensor 
defined by Eq. O. 

Equation (|F2I) possesses an analytic solution in the limit of strong Coulomb interaction e 2 m 2 wo ^S> 1 when the 
ground state wave function reduces to the following asymptotic form 



^o(xi,x 2 ) 



1 



(£l+a=o) 



9?2 



2/2 



(F3) 



22 



where l x = (2H/\/3muJo)^ , and xq — (2e 2 /mwjjs is the 
classical distance between particles, i.e., the distance that 
minimizes the classical energy of two charged particles in 
the ID harmonic potential. The corresponding ground 
state density takes the form of two well separated "blobs" 



n (x) 



1 



A, 



it — 



IT — 



(F4) 



where the size Aoo of the density blobs located at x = 
±xq/2 is 



\ 2y/Smu>o 



(F5) 



The kinetic stress tensor Eq. (1591) for the ground state 
wave function (|F3[) becomes simply proportional to the 
density: 



T (x) 



y/3 + 1 



w n (x) . 



(F6) 



Another technical observation which simplifies calcu- 
lations in the strong interaction limit is that the cross- 
products of the two exponentials in the square brack- 
ets in (|F3j) are irrelevant for the expressions of the type 
^o(x, X2)^o{x' , X2) in the limit of xo^/tocjo 3> 1- For the 
two-particle density matrix entering the nonlocal term in 
(IF2[) this implies the following result 



2* 2 (x,x') = 



t 2 



(F7) 



r 



Simplification of the integral term in the equation of 
motion (|F2I) comes from the fact that in the limit of 
XQ^/muiQ 3> 1 the pair correlation function (|F7[) is peaked 
at \x— x'| ~ x (this keeps the particles at a distance close 
to the classical value). Hence in the integral kernel the 
interaction factor can be approximated as 



d^.w(x — x 1 ) 



2e 2 



2e 2 



Substituting (|F6J) and (JF8J into the (JF2J) and using the 
obvious identity 



n (x) = / dx'2^l(x,x') . 



p'|3 



(F8) we reduce the equation of motion to the following form 



to [w 2 — 2wq] n (x)u(x) — —mui 2 j dx'2^(x, x')u(x') + - — <9 2 [no(x)<9 2 u(x)] — 3a> 



V3 + 1 



d x [n (x)d x u{x)}. (F9) 



r 



Now the following observations are in order: 

(i) The intcgro-differential operator in (|F9I) (in fact in 
the original Eq. (|F2|) ) is symmetric under inversion of x. 
Hence the solutions can be classified by parity w ± (x) = 
±u ± (— x). Therefore it is sufficient to consider Eq. (|F9[) 
only in the region of positive x; 

(ii) In (|F9j) for x > all local terms contain tiq(x) 
which is a narrow Gaussian located at x ~ xq/2. There- 
fore these terms are nonzero only around xq/2; 

(iii) The integral kernel ^(x, x 1 ) for x > is a product 
of two Gaussian peaks, one at x + x' ~ 0, and another at 
x — a;' ~ xo, which confines x to the region of the right 
density blob, x ~ xo/2, and x' to the region of the left 
blob, x 1 ~ — xo/2. 

Therefore for positive x all terms in (|F9| are nonzero 
only in the region x ~ xq/2, while the integration re- 
gion in the nonlocal (interaction) term is confined by the 
Gaussian factors to x' ~ — xq /2. Note that for this reason 



the integral term will contribute with opposite signs to 
the equations of motion for the modes of opposite parity. 

To further simplify the eigenvalue problem in the 
strong coupling limit we proceed as follows, (i) Con- 
sidering Eq. (|F9|) in the region x > we make a shift of 
coordinates x — > x + xo/2, and x' — > x' — xo/2. After 
that because of the Gaussian factors the integration can 
be extended to the whole axis. This completely elimi- 
nates Xo (i.e. the coupling constant) from the problem, 
as it should be in the strong coupling limit; (ii) Go to di- 
mensionless coordinates £ = x/Aqo, and £' = x'/X^; (iii) 
Divide everything by the ground state density no (which 
is simply a Gaussian located at the origin after the above 
shift, no(0 = e-«70O- 

The result of these three steps is the following dimen- 
sionless equation of motion 
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r, iy OO 

W / (v^+l) 2 2 L J ^2^F3 1 /4y ^ vw v ; 

— OO 

To solve the eigenvalue problem fjFlpp we employ the following identities for Hermite polynomials 4 ^ 

e x2 d x (eT x2 d x H k } = -2kH k , e x " d 2 x (e^ d 2 x H k ^ = 4fc(fc - l)H k , (Fll) 

OO , 

1+1 J Me-^V+SM'HuW = (-D fe (^) ' * (*) • (F12) 



From the identities (|Fll[) and (IF12|) we see that Hermite polymomials are the eigenfvmctions of each of three terms in 
the integro-differential operator (both for odd and for even modes) on the right hand side in (|F10|) . The corresponding 
eigenvalues (one should apply the identity (|F12|) with a = are 



"aj (V3 + 1) 2 \V3 + lJ 

where k = 0, 1, 2, . . . is the quantum number labeling the eigenmodes (for every n there are two modes of opposite 
parity). Note that the second, third, and fourth terms on the right hand side of Eq. (|F13|) are, respectively, the 
eigenvalues of the first, second, and the third terms on the right hand side in Eq. (|F10[) . With a little algebra we 
simplify the eigenvalues of Eq. (IF13[) as follows 



1 1/2 

2 + 3v / 3fc + 6fc(/c-l)(2-v / 3)T(-l) fe (2" V3) k ■ (F14) 



In the physical units of length the eigenfunctions in the whole space take the form 

Finally, the displacement eigenmodes normalized by the condition 

dxn {x)u^ {x)uf (x) = Ski (F16) 



can be written as follows 



1 e ^# fc (^^) ± (-l) fc e ( 



£+£o/2 

" fe = Tf^fcT (*-*q/2)2 (x + xq/2)2 ■ ,F1,1 

e A °= + e A ^ 



I 

Equations (|F14[) and (|F17|) give the asymptotically exact of strong correlations. In Sec. IV C we have used this 
solutions of the elastic eigenvalue problem in the limit solutions to control the accuracy of our numerical results. 
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